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Generalized rules for building and flipping clusters in the 
quantum Monte Carlo loop algorithm are presented for the 
XXZ-model in a uniform magnetic field along the Z-axis. As is 
demonstrated for the Heisenberg antiferromagnet it is possible 
from these rules to select a new algorithm which performs 
significantly better than the standard loop algorithm in strong 
magnetic fields at low temperatures. 



I. INTRODUCTION 

The invention of the Loop algorithm JjJ was a major 
breakthrough for Monte Carlo simulations of quantum 
spin systems. This algorithm which is a quantum version 
of the Swendsen-Wang cluster algorithm |2| has many 
desirable features which makes it possible to study large 
systems at low temperatures ||. Among them is that 
the algorithm can be formulated directly in continuous 
imaginary time 0], thus avoiding the need to extrapo- 
late data obtained at finite imaginary time discretization. 
Another is that updates are done in an extended config- 
uration space of spins and new entities called loops. This 
makes big changes in the spin configuration possible in 
a single Monte Carlo step, resulting in very small au- 
tocorrelation times. Furthermore the nonlocal updating 
procedure allows all topological sectors to be sampled. 
For a quantum spin system this means in particular that 
sectors with different magnetization are sampled. This 
is in contrast to most other algorithms which operate at 
fixed magnetization. 

Although excellent for a wide class of models, the Loop 
algorithm does not do well when the Hamiltonian is made 
asymmetric by a uniform magnetic field or a chemical po- 
tential. In these cases autocorrelation times become very 
long at low temperatures and the performance of the al- 
gorithm is lowered drastically. Here it is shown how this 
can be overcome by generalizing the loop algorithm. The 
generalization is obtained by relaxing the condition of 
non-interacting loops, and by taking the magnetic field 
into account in the loop building process. To be spe- 
cific, we consider the nearest neighbor XXZ-model on a 
bipartite lattice in a magnetic field along the Z-axis 
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els as well, such as lattice fermions in the presence of a 
chemical potential. In the next section the generalized 
loop algorithm is presented. Then it is explained how an 
algorithm that performs well in a magnetic field can be 
chosen. The usefulness of this algorithm is demonstrated 
by measuring magnetization curves for the Heisenberg 
antiferromagnet on a dimer, a chain, and on a plane. 
Finally it is shown that the algorithm can also be used 
to determine the critical temperature for the Kosterlitz- 
Thouless transition occuring at finite magnetic fields in 
the Heisenberg antiferromagnet. 



II. THE ALGORITHM 

To explain the algorithm we begin by formulating the 
d dimensional quantum system as a classical system in 
d+1 dimensions. This is done in the standard way |l]] of 
dividing the Hamiltonian into sums of commuting pieces. 
Then a Trotter-Suzuki breakup is performed, and com- 
plete sets of states, which are labelled by their eigenval- 
ues for 5f , are inserted at each time-slice between each 
sum of commuting pieces. The matrix elements are easily 
evaluated and corresponds to interactions around shaded 
plaquettes in a generalized checkerboard pattern. 

As is shown in fig. |l| there are six allowed spin configu- 
rations around a shaded plaquette for the XXZ-model in 
a magnetic field. Other configurations have zero weights 
as for those S z is not conserved along the time direction. 




FIG. 1. The different plaquettes for the XXZ-model in a 
magnetic field. The vertical is the imaginary time direction. 

Because the loop algorithm can be formulated in con- 
tinuous imaginary time it is sufficient to consider the 
limit where the imaginary time spacing At goes to zero. 
In this limit the plaquette weights are 

W (a + ) = !-(£-!) Ar, 
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Despite this choice of model it is expected that the proce- 
dure employed here should apply to other quantum mod- 
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where z is the lattice coordination number. 

The loop algorithm consists of two main steps. The 
first is to build loops. Loop building is a probabilistic 
process where each shaded plaquette p is broken up into 
loop segments G p with a probability P(s p — > s p , G p ), de- 
pendent on the spin configuration s p . Each loop segment 
connects two or four spins. The different types of loop 
segments are shown in fig. |2| 





FIG. 2. The different breakups G into loop segments 
around a shaded plaquette. 

When this is done for all shaded plaquettes, the entire 
space-time lattice will be filled with loops. The second 
step is to flip spins along one or more loops. Because of 
the way the break-ups are constructed, this always re- 
sults in an allowed spin configuration provided all the 
spins around a loop are flipped. The process of flipping 
the spins along a loop is also probabilistic. It is governed 
by the probability Pg{s — > s/) for changing spin con- 
figurations given a particular break-up G for the whole 
lattice. After this second step a new spin configuration is 
generated and one does the measurements and start over 
again. For the whole procedure to satisfy detailed bal- 
ance it is sufficient || that the probabilities are chosen 
such that 
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P G (s -» s/)]Jw(s p ,G p ) = P G (s/~> s)Y[w(s p /,G p ), (4) 

p p 

where G and s are the full loop and spin configuration, 
pieced together by the loop segments G p and plaquette 
spins s p respectively. w(s p , G p ) is the plaquette weight of 
plaquette p in the extended configuration space of both 
spins and loops. The weights w(s p , G p ) must be positive 
definite and satisfy 



^2w(s p ,G p ) = w(t 
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The different loop algorithms described here correspond 
to different choices of these weights. Writing out Eq. (||) 
explicitly we find 

w(a + ) = w(a + , G||) + w(a + ,G x ) + w(a+, G 8 ), (6) 

w(a-) = w(a-,G\\) + w(a_, G x ) + io(o_, G®), (7) 

w(b) = w{b, G=) +w(b, G X )+ W (6,G 8 ), (8) 

w(c) = w(c,Gu) + w(c,G=)+w(c,G®). (9) 



We have set the weights w(a+,G = ), iu(a_, G=), w(b, G||) 
and w(c,G x ) to zero as flipping the spins along one 
loop segment for such configurations leads to a configu- 
ration with zero weight. It is therefore clear that we have 
eight parameters at our disposal. Let us parametrize the 
weights in the following way 
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The remaining four weights are given by Eqs. (||)- 
(|J). Note that although this is a convenient way of 
parametrizing the weights it is not the most general one. 
In selecting the parametrization above we have chosen 
which weights are of order At or 1. 

With this parametrization it is easy to obtain the loop 
building probabilities P(s p — > s p ,G p ) from Eqs. ([}]) and 
(||). To have non- negative weights, all the parameters 
must be greater than or equal to zero and (u+g) < 1^1/2. 

To satisfy detailed balance in the ex- 
tended configuration space Eq. we need the ratios 
w(s p /,Gp)/w(s p ,G p ) which are 
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Given these ratios the flipping probabilities can be gotten 
from 
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III. PARAMETER CHOICES 

There are many possibilities for the choice of param- 
eters, but not all of them lead to efficient ergodic algo- 
rithms. To minimize autocorrelation times one must in 
particular ensure that the loops generated have a reason- 
able chance of being flipped. This means that the ratios 
in Eqs. (Es|)-(E^) should be as close to unity as possible. 

Let us first consider H = 0. The standard loop algo- 
rithm is constructed such that all the ratios Eqs. (|l8|)- 
(pq) are one, and by minimizing the weights w(x, G®): 



u = v = e = f = g = h = 0, s = t= -y-, 
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In particular the nonzero parameters for the Heisenberg 
antifcrromagnet corresponds to, uq = vo = J/2, and for 
the XY-model sq — to = u n = v n = J/4. For Ising 
anisotropy, \J X \ < |J Z |, certain G® break-ups must be 
included as otherwise some weights will be negative. For 
extreme anisotropy \J X \ = the model is the classical 
Ising model and the world- lines are all straight, sq = to = 
vq = 0. In this limit the standard loop algorithm above 
is the Swendsen-Wang algorithm for the Ising model. 

Now consider H ^ 0. In this case it is not possible 
to set all of the ratios Eqs. (p^|)-(p5|) to unity for any pa- 
rameter choices. With the parameter choices for the stan- 
dard loop algorithm these ratios are only unity for loops 
which do not change the magnetization when flipped. For 
loops that can change the magnetization the total ratio 
of weights is exp(— pHAM), where AM is the change 
in magnetization caused by flipping the loop. This leads 
to autocorrelation times that increase exponentially with 
f3H. The reason is that the magnetic field is not taken 
into account in the loop building process. The process 
of changing the magnetization is a competition between 
loosing Zeeman energy and gaining exchange energy. As 
the exchange energy is gained in the loop building pro- 
cess, it is inefficient to build the loops as if the magnetic 
field was absent. What happens in the standard loop al- 
gorithm is that the number of loops generated which can 
change the magnetization is very small. 

An interesting observation is that one can construct an 
algorithm where only loops which can change the mag- 
netization are generated. This choice is 



which means that the only break-ups allowed are of the 
diagonal type. This algorithm is ergodic, and in the clas- 
sical Ising limit, \J X \ = 0, it corresponds to the standard 
local Ising model algorithm in a magnetic field. One can 
now ask for which magnetic field this algorithm mini- 
mizes the autocorrelation times. Eq. (18) is unity for 
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The important observation is that this value of H is the 
saturation field, where almost all plaquettes are of type 
a + . At this field it is therefore not important for the 
performance of the algorithm that Eqs. ( |l9| ) and j20| ) 
deviate from one. 

It is then natural to choose an algorithm valid for all 
H which interpolates between the standard algorithm at 
H = and the above at the saturation field. For the 
Heisenberg antiferromagnet in a magnetic field we thus 
propose the following algorithm 
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For H > Jz we use the same algorithm as for H = Jz. 
Eq. @ implies that Eqs. becomes unity, 

whereas Eqs. ©-(pOl) do not. 



IV. NUMERICAL RESULTS 

The algorithm was first tested by measuring magne- 
tization of a dimer or two-site 5=1/2 antiferromag- 
netic chain. As shown by Kashurnikov et al. Q the inte- 
grated autocorrelation time for the standard loop algo- 
rithm increases exponentially with (3H . We have verified 
this measuring the integrated autocorrelation time as de- 
scribed in Ref.l, appendix B. For the new algorithm the 
dimer integrated autocorrelation time is very small (< 2) 
down to the lowest temperatures measured (T=.005J) 
and the magnetization measured agrees excellently with 
exact results. 

Fig. H shows the magnetization per spin of a 64-site 
antifcrromagnetic Heisenberg chain, and illustrates the 
improvement over the standard Loop algorithm. The re- 
sults of the modified and standard loop algorithm were 
obtained using the same number of equilibration and 
measurement steps (10 6 ), and the lines are exact results 
obtained using Bethe Ansatz 0. It is clear that, in con- 
trast to the modified algorithm, results obtained using 
the standard loop algorithm have not converged for high 
magnetic fields. Close inspection of the data at the lowest 
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temperature reveals that the results of the modified algo- 
rithm deviate slightly from the exact results at interme- 
diate magnetic fields. This deviation which is statistical 
is caused by increased autocorrelation times which arises 
because Eqs. (|l9|)-(p0|) deviate from unity concomitant 
with the presence of a significant fraction of a_ plaque- 
ttes at these fields. For the 64-site chain we measured the 
integrated autocorrelation time to be a maximum 3 • 10 
steps at H/J — 1.3, going down to about 60 steps at low 
and high fields. It is quite conceivable that a different 



interpolation scheme than the one chosen in Eq. (37) can 



reduce these autocorrelation times at intermediate fields. 
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FIG. 3. Magnetization per spin for a spin chain with 64 
sites. 

Fig. ^ shows the full magnetization curve of a 64x64 
square lattice antiferromagnet. Typical runs involved 
10 7 steps for equilibration and measurements. The inset 
shows the high field behavior for very low temperatures 
on the same lattice. The statistical errors, taking into 
account the autocorrelation times, are smaller than the 
symbol size. For the lowest temperature the integrated 
autocorrelation times reached a maximum of 4 • 10 4 steps 
at H/J=2.4 going down to about 4 • 10 3 steps at low and 
high magnetic fields. 



The Heisenberg antiferromagnet in a magnetic field un- 
dergoes a Kosterlitz-Thouless transition at finite temper- 
atures. The transition temperature has previously been 
obtained for weak mag netic fields; H < .2 J §. Fig. | 
shows the helicity modulus T, which is the normalized 
free energy change due to phase twists in the x-y-plane, 
and which is proportional to the squared spatial winding 
number of world- lines Q , as a function of temperature 
for four different system sizes at H = 3.95 J. From a fi- 
nite size analysis [ [To| the transition temperature is found 
to be T c = .020(5) J. Here a single-cluster [[TJ implemen- 
tation of the algorithm is used. It is expected that more 
precise estimates for T c can be obtained using a multi- 
cluster implementation. 




FIG. 5. Helicity modulus as function of temperature 
for two different system sizes. The line is the Koster- 
litz-Thouless-Nelson critical line. 

The author wishes to thank P.A. Lee for useful discus- 
sions and V. Chudnovsky and U.-J. Wiese for providing 
source code for the continuous time loop algorithm. 
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FIG. 4. Magnetization per spin for a plane with 64 x 64 
sites. 
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